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ABSTRACT 

The structural and dynamical properties of star clusters are generally derived by 
means of the comparison between steady-state analytic models and the available ob¬ 
servables. With the aim of studying the biases of this approach, we fitted different 
analytic models to simulated observations obtained from a suite of direct N-body sim¬ 
ulations of star clusters in different stages of their evolution and under different levels 
of tidal stress to derive mass, mass function and degree of anisotropy. We find that 
masses can be under/over-estimated up to 50% depending on the degree of relaxation 
reached by the cluster, the available range of observed masses and distances of radial 
velocity measures from the cluster center and the strength of the tidal field. The mass 
function slope appears to be better constrainable and less sensitive to model inad¬ 
equacies unless strongly dynamically evolved clusters and a non-optimal location of 
the measured luminosity function are considered. The degree and the characteristics 
of the anisotropy developed in the N-body simulations are not adequately reproduced 
by popular analytic models and can be detected only if accurate proper motions are 
available. We show how to reduce the uncertainties in the mass, mass-function and 
anisotropy estimation and provide predictions for the improvements expected when 
Gaia proper motions will be available in the near future. 

Key words: methods: numerical - methods: statistical - stars: kinematics and dy¬ 
namics - globular clusters: general 


1 INTRODUCTION 

The internal kinematics of star clusters offers a wealth of 
information on their present day dynamical status and pro¬ 
vides precious traces of the past history of evolution and 
interaction with the host galaxy of these stellar systems. 
In old and relatively low-mass stellar systems (like open 
and globular clusters) the half-mass relaxation time is of¬ 
ten shorter that their ages and processes like kinetic energy 
equipartition, mass segregation and core collapse can be at 
work and leave signatures in the phase-space distribution of 
their stars. The understanding of these processes is crucial 
to properly model the dynamics of these stellar systems and 
to derive their global characteristics (total mass, mass func¬ 
tion, degree of anisotropy, etc.) from observations performed 
in spatially restricted region of the cluster. 

The easiest way to model the internal kinematics of star 
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clusters is through analytic models. Alternative approaches 
include Jeans modelling (Ogorodnikov, Nezhinskii & Os- 
ipkov 1978; Mamon & Boue 2010), Schwarzschild’s orbit 
superposition method (Schwarzschild 1979; Vasiliev 2013) 
and comparison with Monte Carlo and N-body simulations 
(Giersz et al. 2008, 2009, 2011; Heggie et al. 2014). However, 
the former two approaches do not ensure physically mean¬ 
ingful solutions (i.e. the corresponding distribution function 
could be negative), while the latter two require an enormous 
computational effort and have started to become feasible 
for modelling rich star clusters only in recent years (e.g., see 
Heggie et al. 2014). Analytic models are generally defined by 
distribution functions depending on constants of the motion, 
and assume the cluster is in a steady-state and in equilib¬ 
rium with the surrounding tidal field. Because star clusters 
are good approximation of collisions-dominated systems we 
have a relatively advanced understanding of the distribu¬ 
tion of their stars in phase-space from theory and numerical 
simulations, and the choice of distribution function-based 
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models is justified. The most popular model of this kind is 
the King (1966) model which proved to be quite effective in 
reproducing the surface brightness profiles of many globular 
clusters (GCs), open clusters and dwarf galaxies (Djorgovski 
1993; McLaughlin & van der Marel 2005; Carballo-Bello et 
al. 2012; Miocchi et al. 2013). A generalization of this model 
accounting for radial anisotropy and a degree of equiparti- 
tion among an arbitrary number of mass components has 
been provided by Gunn & Griffin (1979; see also Da Costa 
& Freeman 1976; Merritt 1981). The underlying assumptions 
of these models (i.e. the functional dependence of the distri¬ 
bution function on the integrals of motion and masses), al¬ 
though relying on a physical basis, are only arbitrary guesses 
to model the result of the complex interplay among many 
physical processes. Unfortunately, the validity of these as¬ 
sumptions is not easy to be tested with real data since it 
would require the determination of complete radial profiles, 
radial velocities and relative proper motions of large sam¬ 
ples of stars with different masses. However, ground based 
photometric observations are not able to resolve faint low- 
mass stars in the cores of the most concentrated globular 
clusters (GCs), and Hubble Space Telescope (HST) obser¬ 
vations are often available only for a limited portion of the 
cluster extent. Therefore, a single density profile is often cal¬ 
culated by mixing surface brightness estimates (dominated 
by the contribution of red giant stars) in the innermost re¬ 
gion with main sequence (MS) star counts in the outskirts. 
For the same observational limits, information on the shape 
of the low-mass end of the mass function (MF) are avail¬ 
able only for a few nearby GCs. Moreover, accurate line-of- 
sight (LOS) velocities can be derived through spectroscopic 
analyses only for the most massive red giant branch (RGB) 
stars. An even more complex task is the derivation of proper 
motions: to date, the typical uncertainty of the available as¬ 
trometric surveys and their limiting magnitude allowed the 
estimate of relative proper motions only for the brightest 
stars of the closest GCs (McLaughlin et al. 2006; Ander¬ 
son & van der Marel 2010). These data, while providing an 
indication of the systemic cluster motion, are not accurate 
enough to be used for studies on their internal dynamics. 
Only very recently proper motions with the right accuracy 
derived from extensive HST surveys are becoming available 
(Bellini et al. 2014; Watkins et al. 2015). In the absence of 
suitable observational data, Trenti & van der Marel (2013) 
studied the effect of two-body relaxation on the projected 
velocity dispersion of groups of stars with different masses in 
a set of N-body simulations. They noted that the simulated 
clusters never reach complete kinetic energy equipartition 
(a result already found previously by Giersz & Heggie 1997 
and confirmed on real GCs by preliminary results by Bellini 
et al. 2013) and argued that the widely used King-Michie 
models could be inadequate to model these stellar systems. 

The possible inadequacy of these models could have se¬ 
rious consequences on the estimate of many intrinsic prop¬ 
erties of GCs obtained so far. For example, since a GC is 
a self-gravitating system in virial equilibrium, its velocity 
dispersion can be used as a proxy for its mass. For this 
reason, many studies have been devoted in the past to the 
determination of dynamical masses by means of the compar¬ 
ison of the available photometric and kinematical observ¬ 
ables (mainly projected density/surface brightness profiles 
and line-of-sight radial velocities of RGB stars) with ana¬ 


lytic models (Mandushev, Staneva & Spasova 1991; Pryor 
& Meylan 1993; Lane et al. 2010; Zocchi, Bertin & Varri 
2012; Kimmig et al. 2014). Because of the above mentioned 
lack of information on the radial distribution and kinemat¬ 
ics of low-mass stars, these studies estimated masses through 
the comparison of the available data with single-mass mod¬ 
els (such as King 1966 models) using RGB stars as tracers 
of the potential. This approximation, however, is incorrect 
and can lead to significant biases in those clusters where a 
large number of collisions allowed an effective exchange of 
kinetic energy among stars with different masses (Shana¬ 
han & Gieles 2015, hereafter SG15). Indeed, massive stars 
tend to transfer kinetic energy in collisions with less mas¬ 
sive ones, becoming kinematically cold and moving on less 
energetic orbits at small radii. So, in these clusters both the 
spatial distribution and the velocity dispersion of RGB stars 
significantly differ from those of the other less massive stars 
which contribute to the largest fraction of the cluster mass. 
The situation is further complicated by the presence of dark 
remnants (white dwarfs, neutron stars and black holes) and 
binary stars contributing to the cluster potential but whose 
fraction is hard to be estimated. Because of these biases, 
the comparison of the derived dynamical masses with those 
estimated converting the cluster integrated luminosity into 
mass (often referred as ’’luminous mass”) led to puzzling 
and conflicting results. Indeed, the first studies indicated 
that dynamical masses were ~25% smaller than luminous 
ones (McLaughlin & van der Marel 2005; Strader et al. 2009, 
2011). The above discrepancy is partly due to the adoption 
of RGB stars as dynamical tracers of the cluster potential, 
and to the neglecting of dynamical processes like the prefer¬ 
ential losses of low-mass stars altering the mass-to-light ra¬ 
tio used to derive the luminous masses (Kruijssen & Mieske 
2009). The above effects have been accounted for in a recent 
analysis by Sollima, Bellazzini & Lee (2012) who compared 
deep HST photometric data and an extensive survey of ra¬ 
dial velocities for six Galactic GCs with a set of multi-mass 
King-Michie models, simultaneously fitting the luminosity 
function estimated in the innermost cluster region, the sur¬ 
face brightness and the velocity dispersion profile of RGB 
stars to derive global masses and MFs. Adopting the recipe 
of mass segregation of these models, they found that dy¬ 
namical masses are on average ~40% larger than luminous 
ones i.e. in the opposite direction of what found in previous 
analyses. In this case, while a large portion of this discrep¬ 
ancy can be due to an uncertain assumption of the fraction 
of dark remnants, a significant contribution could be due to 
the inadequacy of the adopted models in reproducing the 
degree of mass segregation of the analysed clusters. 

Another aspect connected to the modelling of mass seg¬ 
regation of GCs regards the derivation of their global MF. 
Indeed, thanks to the proximity of many GCs and the high 
angular resolution of HST it has been possible to sample 
their luminosity function down to the hydrogen burning 
limit (Piotto, Cool & King 1997; Piotto & Zoccali 1999; 
Paust et al. 2010). These kind of observations are however 
limited to a restricted portion of the cluster, generally close 
to the half-mass radius (where the effects of mass segrega¬ 
tion are minimized; Pryor, Smith & McClure 1986; Vesperini 
& Heggie 1997; Baumgardt & Makino 2003) or the cluster 
center (see e.g. the ”ACS globular clusters treasury project ”; 
Sarajedini et al. 2007). Because of mass segregation, to de- 


© 2015 RAS, MNRAS OOO.HlfHl 


Biases in GCs parameters estimation 3 



t [Gyr] t [Gyr] 


Figure 1. Lagrangian radii evolution of the N-body simulation W5rhll.5R8.5 (left panel) and W5rhlR8.5 (right panel). The radii 
containing the 1%, 2%, 5%, 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 90% of the cluster mass are marked. The time at which selected 
snapshots are extracted are marked with dashed lines. 


rive the global MF from the local one a correction based on 
multi-mass modelling is necessary. This task has been per¬ 
formed by Paust et al. (2010) and Sollima et al. (2012) who 
found that the MFs of their samples of GCs can be well re¬ 
produced by power-laws with slopes varying in a wide range 
(—1.7 < clmf < —0.2). Any possible bias in the adopted 
prescription for the distribution of masses would translate 
into a systematic error in the estimated global MF. 

Finally, it is not clear if the velocity distribution of GCs 
stars is isotropic or if it presents some radial/tangential bias. 
Physical reasons at the origin of such anisotropy can be 
found in their mechanism of formation (possibly related to 
the collapse of their original gas cloud; Lynden-Bell 1967), to 
the collisions occurring in their central regions (Lynden-Bell 
& Wood 1968; Spitzer & Shull 1975) or to the interaction 
with the Milky Way (Oh & Lin 1992). The determination 
of anisotropy in GCs is extremely complex since it requires 
the knowledge of transverse velocities for a large sample of 
stars to allow an unambiguous detection, and detailed stud¬ 
ies have been possible only very recently for the innermost 
regions of nearby GCs (Watkins et al. 2015). King & Ander¬ 
son (2001) suggested that the dependence of the degree of 
anisotropy on mass could be different from what is predicted 
by King-Michie models. 

The advent of the Global Astrometric Interferometer 
for Astrophysics (Gaia) will change the picture by provid¬ 
ing accurate distances and proper motions for hundreds of 
stars in a large sample of nearby GCs (Pancino, Bellazzini 
& Marinoni 2013). Indeed, the unprecedented astrometric 
accuracy of Gaia will provide proper motions with an accu¬ 
racy comparable to that of LOS radial velocities obtainable 
through high resolution spectroscopy. These measures, com¬ 
plemented by those soon available from the HST analyses 
in the central regions of GCs, could open the way for a new 
golden age for studies on the internal dynamics of GCs. 

In this paper we will compare the most widely used 
analytic models with a suite of N-body simulations apply¬ 
ing the most widely used technique to determine mass, MF 


and anisotropy with the aim of studying the biases in the 
determination of these quantities linked to the ability of an¬ 
alytic models to reproduce the correct phase-space struc¬ 
ture of the simulated stellar systems. In Sect. 2 the adopted 
analytic models and the performed N-body simulations are 
described. In Sect. 3 we describe the method adopted to de¬ 
rive the best-fit masses and MFs. In Sect. 4 the results of 
the comparison between the true and the estimated masses 
and MFs are presented. Sect. 5 is devoted to the simulation 
of Gaia observations of GCs including realistic errors and 
radial sampling efficiency to test the feasibility of studies on 
GCs anisotropy. We finally summarize our results in Sect. 6 


2 SIMULATIONS AND MODELS 
2.1 N-body simulations 

The N-body simulations considered here have been per¬ 
formed using the collisional N-body codes NBODY4 and 
NBODY6 (Aarseth 1999) and are part of the surveys pre¬ 
sented by Baumgardt & Makino (2003) and Lamers, Baum- 
gardt & Gieles (2013)0- Each simulation contains 131072 
particles with no primordial binaries, although a small num¬ 
ber of binaries and triples form through tidal capture during 
the cluster evolution. In this configuration, the total clus¬ 
ter mass is 71236.4 Mg. Particles were initially distributed 
following a King (1966) model with central dimensionless 
potential Wo = 5, regardless of their masses. Two simu¬ 
lations with different half-mass radii have been run (with 
r h = 1 and 11.5 pc, hereafter referred to as W5rhlR8.5 and 
W5rhll.5R8.5, respectively). Particle masses taken from a 
Kroupa (2001) MF with a lower mass limits of 0.1 Mg 
and an upper mass limit of 15 and 100 Mg, for the 
W5rhll.5R8.5 and W5rhlR8.5 simulation respectively. The 

1 both simulations are available at the ” Gaia Challenge” webpage 
http://astrowiki.ph.surrey.ac.uk/dokuwiki 
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cluster moves within a logarithmic potential having circu¬ 
lar velocity iWc = 220 km/s, on a circular orbit at a dis¬ 
tance of 8.5 kpc from the galactic center. The correspond¬ 
ing initial Jacobi radius is rj = 61.15 pc, i.e. equal to the 
one of the W5rhll.5R8.5 simulation. Because of their dif¬ 
ferent Roche lobe filling factors, the tidal field affects the 
two simulations in extremely different ways. Moreover, the 
initial half-mass relaxation time is significantly longer in 
model W5rhll.5R8.5 ( t r h = 4.97 Gyr) with respect to model 
W5rhlR8.5 ( t r h = 0.12 Gyr). So, while the mass of the con¬ 
sidered simulations are significantly smaller than those of 
real GCs, they bracket the majority of globular clusters in 
terms of both their half-mass relaxation time and relative 
strength of the tidal field ( rj/m ). Stellar evolution has been 
modeled using the fitting formula of Hurley, Pols & Tout 
(2000) assuming a metallicity Z=0.001. Mass lost during 
stellar evolution is assumed to be immediately lost from the 
cluster. No kick velocity has been added to stellar remnants 
at their birth. A certain retention fraction of newly formed 
black holes and neutron stars has been assumed: in simula¬ 
tion W5rhll.5R8.5 no black holes are present and all neu¬ 
tron stars are assumed to be retained at their birth, while in 
simulation W5rhlR8.5 we assumed that 10% neutron stars 
and black holes are retained. The temperature and luminos¬ 
ity of each star is recorded at each time-step according to its 
evolutionary stage. Simulations have been run until cluster 
dissolution, defined to be the time when 95% of the initial 
mass has been lost. The evolution of the Lagrangian radii 
of the two simulations is shown in Fig. |T] After the initial 
expansion driven by stellar evolution mass-loss, simulation 
W5rhlR8.5 quickly undergoes core-collapse (after ~5 Gyr), 
and continues its evolution expanding slowly. Instead, be¬ 
cause of its longer relaxation time, simulation W5rhll.5R8.5 
follows a slow evolution losing a large fraction of its stars be¬ 
cause of the strong tidal field. 


2.2 Analytic models 


Selected snapshots of the above N-body simulations have 
been extracted and their projected properties (density and 
LOS velocity dispersion profiles) have been compared with 
a set of King-Michie models (Gunn & Griffin 1979). These 
models are constructed from a lowered-Maxwellian distribu¬ 
tion function made by the contributions of H mass groups 
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where E and L are, respectively, the energy and angular 
momentum per unit mass, v r and vt are the radial and tan¬ 
gential components of the velocity, the effective potential xp 
is the difference between the cluster potential <f> at a given ra¬ 
dius r and the potential at the cluster tidal radius xp = <p—<pt., 
Ai and ki are scale factors for each mass group, and cjk is a 
normalization term. Note that the above distribution func¬ 
tion allows for various levels of radial anisotropy but does 


not allow tangential anisotropy. The number density and the 
radial and tangential components of the velocity dispersion 
of each mass group are obtained by integrating Eq. |Tj 
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Equations [2] can be written in terms of dimensionless 
quantities by substituting 
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is the core radius (King 1966). The potential at each radius 
is determined by the Poisson equation 


SJ %p = 4nGp 


(3) 


Equations [5] and 0 have been integrated after assuming as 
a boundary condition a value of the dimensionless potential 
at the center Wo outward till the radius ry at which both 
density and potential vanish. The total mass of the cluster 
M is then given by the sum of the masses of all the groups 
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As a last step, the above profiles have been projected onto 
the plane of the sky to obtain the surface density 
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the velocity dispersion along the line-of-sight, and in the 
projected radial and tangential directions in the plane of 
the sky 
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The MF of unevolved stars in an annulus defined at a given 
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distance R and with a width A R is given by 
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where p,i = 1 — N[ ernrl /Ni is the fraction of unevolved stars 
and NT ernrl j s the number of remnants in the i-th mass bin. 

In the above models, for any choice of the Ai co¬ 
efficients, the shape of the density and velocity disper¬ 
sion profiles are completely determined by the parameters 
( Wo,f a ,Ni ) while their normalization is set by the pair of 
parameters ( r c , M). 

The dependence on mass of the coefficients A, deter¬ 
mines the degree of mass segregation of the cluster. In the 
formulation by Gunn & Griffin (1979) Ai oc rrn which im¬ 
plies that more massive stars are kinematically colder than 
less massive ones. Note however that, according to Eq.[2]and 
assuming isotropy, the actual squared velocity dispersion of 
the i-th mass group is 

2 = 6 a 2 K 7(7/2, AjW) 
ai 5 Ai 7(5/2, AiW) 

where 7 ( 0 , 6 ) = f b x a ~ 1 e~ x dx is the lower incomplete 
gamma function. According to the above equation, if Ai oc 
rrii complete equipartition is achieved only in the limit 
W —> + 00 . So, these models do not assume kinetic energy 
equipartition neither locally nor globally (see also Merritt 
1981). A slight modification to the Gunn & Griffin (1979) 
assumption can be made to allow intermediate levels of re¬ 
laxation, by simply assuming Ai oc mf. In this case a = 0 
models correspond to single-mass King (1966) models, a = 1 
to Gunn & Griffin (1979) models. We note that for a —» +00 
the distribution function of low-mass stars tends asymptot¬ 
ically to f(E) oc —AiE/2 <tk i.e. there is a maximum degree 
of kinetic energy equipartition reachable by these models. 
Moreover, extrapolation at values of a >> 1 produces un¬ 
realistic density and velocity dispersion profiles. 


3 METHOD 

The aim of this paper is to check the ability of analytic mod¬ 
els to reproduce the main kinematical properties of simu¬ 
lated star clusters. For this purpose, we used selected snap¬ 
shots of the simulations described in Sect. ED to produce 
mock observations, by considering only those information 
generally available to observers. In particular, for each anal¬ 
ysed snapshot we extracted: 

• the projected number density profile of the brightest 
stars; 

• the MF of unevolved stars (N° bs (Rd.)) estimated in a 1 
pc-width annulus around a given distance from the cluster 
center (hereafter referred as the deep field range ); 

• the LOS velocities of the brightest stars. 

Projected distances and velocities have been calculated by 
projecting 3D positions and velocities along a random di¬ 
rection. We considered only stars whose projected distance 


lies within the Jacobi radius, to avoid the inclusion of tidal 
tails in our sample. As described in Sect. 12.21 analytic mod¬ 
els are constructed from the superposition of different mass 
bins. Here we considered 8 evenly spaced mass bins rang¬ 
ing from 0.1 Mq to the mass at the RGB tip ( Mu p ), plus 
an additional bin containing all the remnants more mas¬ 
sive than Mtip. The particles of the simulation have been 
divided in the above defined mass bins according to their 
masses. As we are interested in systematic biases (and not 
in random uncertainties) we have not added observational 
errors and incompleteness effects to our data. Moreover, we 
assumed that the fraction (pi) and the mass distribution of 
dark remnants (fV/ emn ) are known a priori by the observer. 
This last assumption is far from reality and in real cases this 
uncertainty can even represent the largest source of system- 
atics (see Sollima et al. 2012). However, as our purpose is 
to determine the impact of the adopted criterion for relax¬ 
ation we neglect this issue. The impact of the uncertainties 
on dark remnants will be addressed in a forthcoming paper 
(Peuten et al., in preparation). 

The best fit to the above defined quantities has been 
performed using a Markov-Chain Monte Carlo (MCMC) al¬ 
gorithm nested into an iterative procedure (see Sollima et 
al. 2012). We performed the following steps: 

(i) We adopted as a first guess of the global MF (Ni) the 
sum of the local MF of unevolved stars measured in the deep 
field range ( N° b3 ) and that of remnants (N( ernn )-, 

(ii) the MCMC algorithm samples the parameter space 
to find the pair of parameters (Wo, r c ) minimizing the 
Kolmogorov-Smirnov penalty function i.e. the maximum ab¬ 
solute difference between the observed and the predicted 
normalized cumulative distribution of particles in the con¬ 
sidered mass bins (excluding remnants; see Eq.0. The prob¬ 
ability associated to the Kolmogorov-Smirnov test Pks is 
recorded; 

(iii) The MF predicted by the model in the deep field 
range (Nfi lod ) is calculated using Eq. [5] and a new guess 
(N() of the global MF of unevolved stars is made using the 
relation 


N' = Ni 


( N? s y 

^ jymod J 


where the exponent 0.5 is a damping factor to avoid diver¬ 
gence. 

Steps (ii) and (iii) are repeated until the maximum 
variation of the estimated MF in two subsequent iterations 
among the various mass bins (maa:|l — N[/Ni |) falls below 
1 %. 

The mass of the model has been finally derived by min¬ 
imizing the log-likelihood function 


1 / _ .2 

InL = IuPks - - V 7 ~ ^ , + ln(a LOS (Ri)) 

J i=1 a LOS\Ei) 


(7) 


Where v % is the LOS radial velocity of the i-th star 
belonging to the most massive group. An additional cycle 
can be added to the above procedure to search for the value 
of the anisotropy radius f a maximizing the above defined 
merit function. 
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Figure 2. Cumultive radial distribution of different mass groups in the selected snapshots of the performed N-body simulations (solid 
lines). The prediction of the best-fit King-Michie model (assuming a = 1) are marked with dashed lines. 


4 RESULTS 

As a first test we selected two snapshots for each simulation 
in selected moments of the cluster evolution and compared 
the radial distribution, velocity dispersion and anisotropy 
of different mass groups with those predicted by the corre¬ 
sponding best-fit models. In particular, for the less evolved 
simulation W5rhll.5R8.5 we selected a snapshot at 1 Gyr 
(after the stellar evolution-driven expansion and before a 
significant level of relaxation has been reached) and a snap¬ 
shot at 5 Gyr (t ~ 0.8 t r h{t)). For simulation W5rhlR8.5 
snapshots at 5 Gyr (close to core collapse) and 11 Gyr (af¬ 
ter several half-mass relaxation times; t ~ 9 t r h(t)) have been 
selected. 

In Fig. [5] the cumulative radial distribution of four dif¬ 
ferent mass groups are compared with those predicted by 
the best-fit King-Michie models assuming a=l. As expected, 
in the early stage of evolution of simulation W5rhll.5R8.5 
different mass groups have still very similar distributions. 
In this case, the best-fit multi-mass model overpredicts the 
level of relaxation predicting low-mass stars distributed on a 
more extended structure than they really are. On the other 
hand, after 5 Gyr the cluster reaches the level of relaxation 
adopted by the model, which reproduces extremely well the 
radial distribution of all the mass groups. A different situ¬ 
ation occurs in simulation W5rhlR8.5: in both the selected 


stages of evolution the actual mass segregation is slightly 
larger than that predicted by King-Michie models. 

The same conclusion can be reached by analysing the 
ratio between the 3D velocity dispersion^] of different mass 
bins (see Fig. [3]). In a relaxed stellar system kinetic energy 
equipartition implies that at each radius the velocity disper¬ 
sion of stars scales with their masses as a oc mT 1 ^ 2 . So, while 
at the beginning of the simulation stars of different masses 
have the same dispersion, as relaxation proceeds massive 
stars tend to became kinematically colder than less massive 
ones and the ratio of their velocity dispersions tends to the 
square root of their mass ratios. This process is faster in 
the center of the cluster where a large number of collisions 
accelerates this process. An inspection of Fig. [3] shows that, 
as already mentioned in Sect. 12.21 both N-body simulations 
and analytic models are still far from complete kinetic en¬ 
ergy equipartition even in the cluster center. Again, models 
overpredict relaxation in the first selected snapshot of simu¬ 
lation W5rhll.5R8.5 while they well reproduce the velocity 
dispersion ratio in all the other considered snapshots of both 
simulations, except for the least massive bin which appears 


2 Hereafter, in the comparison between the velocity dispersion of 
models and N-body snapshots the systemic motion of the cluster 
is subtracted i.e. cr 2 = (v 2 ) — {v ) 2 . 
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Figure 3. Ratio between the velocity dispersion of the most massive and those of less massive bins as a function of distance from 
the cluster center for the selected snapshots of the performed N-body simulations. The prediction of the best-fit King-Michie models 
(assuming a = 1) are shown with solid lines. The values corresponding to kinetic energy equipartition are marked with arrows. 


over-relaxed in simulation W5rhlR8.5 than what predicted 
by analytical models. 

In Fig.[4]the anisotropy parameter /3 = 1 — at /2crj? of the 
various mass groups is shown as a function of the distance 
from the cluster center. A positive value of /3 indicates a bias 
toward radial orbits (radial anisotropy ), while a negative 
value indicates the prevalence of tangential orbits (tangential 
anisotropQ). It is apparent that the two simulations present 
opposite trends: in both the considered snapshots of simu¬ 
lation W5rhll.5R8.5 tangential anisotropy which increases 
toward the outer cluster regions is apparent in all mass bins 
with the same indistinguishable amplitude. Instead, in sim¬ 
ulation W5rhlR8.5 all mass groups quickly develop radial 
anisotropy. After many relaxation times, however, this trend 
tends to be erased in the less-massive stars whose orbits be¬ 
came again isotropic, while a significant radial anisotropy 
persists in the most massive stars. The physical reason of 
this behaviour is likely due to the different tidal stress ex¬ 
perienced by the two simulations: in simulation W5rhlR8.5 
collisions occurring in the cluster center kick stars in the 
cluster halo on radial orbits. This process, already studied 

3 We define tangential anisotropy as a bias toward tangential 
orbits with random directions. Rotation, characterized by ordered 
motions in a single direction, is not analysed here. 


in past studies by Lynden-Bell & Wood 1968 and Spitzer & 
Shull 1975, appears to affect stars of different masses with 
the same efficiency. In the initial stage of this simulation 
the cluster is very compact and almost unaffected by the 
tidal field. As evolution proceeds, the cluster expands and 
the kinematically hottest low-mass stars reach the bound¬ 
ary of the Roche volume. In this case, stars on radial orbits 
reach this boundary with positive velocity and escape more 
efficiently from the cluster. This process produces the ob¬ 
served increasing trend of radial anisotropy with mass. The 
same effect is at work during the entire evolution of simu¬ 
lation W5rhll.5R8.5: in this case, the cluster feels a strong 
tidal field and stars on tangential orbits are preferentially 
retained regardless of their masses. The comparison with 
King-Michie models has been done only for the snapshots of 
the W5rhlR8.5 simulation. Indeed, these models do not ac¬ 
count for tangential anisotropy present in the W5rhll.5R8.5 
simulation. It is apparent that King-Michie models predict 
an increasing anisotropy for more massive stars. So, while 
they reproduce fairly well the snapshot at 11 Gyr, they can¬ 
not provide a satisfactory representation of the cluster in its 
early stages. 

Summarizing, King-Michie models with a = 1 (i.e. 
Ai oc mi as in the formulation by Gunn & Griffin 1979) ap¬ 
pear to qualitatively reproduce the radial distribution and 
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Figure 4. Anisotropy parameter as a function of the distance from the cluster center for various mass groups in the selected snapshots 
of the performed N-body simulations. The prediction of the best-fit King-Michie models (assuming a = 1) are marked with solid lines. 


relaxation status of the performed N-body simulations af¬ 
ter a timescale comparable to the half-mass relaxation time. 
On the other hand, the formalism of these models does not 
allow to reproduce the general characteristics of the veloc¬ 
ity anisotropy in all the stages of evolution of the simulated 
clusters. 

To quantify the ability of these models to reproduce 
the estimated mass and MFs we fitted models to a large 
number of snapshots of the two considered simulations under 
different assumptions on the available range of stars used to 
compute the density profile and the location of the deep field 
used to estimate the local MF. In particular, we considered 
the following fits: 

(i) single-mass King (1966) models (a = 0) fitting the 
number density profile constructed with the 4 most massive 
bins (red lines in Fig.s0and Q: 

(ii) single-mass King (1966) models (a = 0) fitting the 
number density profile constructed with only the most mas¬ 
sive bin (magenta lines); 

(iii) multi-mass King-Michie models (a = 1) fitting the 
number density profile constructed with the 4 most massive 
bins and adopting a deep field range centered around the 
projected half-mass radius (blue lines); 

(iv) multi-mass King-Michie models (a = 1) fitting the 
number density profile constructed with only the most mas¬ 


sive bin and adopting a deep field range centered around the 
projected half-mass radius (cyan lines); 

(v) multi-mass King-Michie models (a = 1) fitting the 
number density profile constructed with the 4 most massive 
bins and adopting a deep field range centered in the cluster 
center (green lines); 

(vi) multi-mass King-Michie models fitting the number 
density profile constructed with the 4 most massive bins 
and a value of a chosen to simultaneously reproduce the MF 
measured in two deep field ranges centered at the projected 
half-mass radius and in the cluster center, respectively (grey 

lines J3- 

Only isotropic models have been considered since i) 
as shown above, King-Michie models do not reproduce the 
qualitative trend of anisotropy with mass, ii) they do not 
include a treatment for tangential anisotropy occurring in 
some stage of evolution, and iii) in absence of proper mo¬ 
tion information the signal of radial anisotropy (character¬ 
ized by an increase of the central LOS velocity dispersion) 
can be mimicked by the presence of massive objects not 
(properly) accounted in the model (e.g. a large number of 


4 Because of the unrealistic behaviour of models with large values 
of a (see [2l2t we performed this comparison only with snapshots 
with a best-fit value of a <1. 
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Figure 5. Ratio between the masses estimated by analytic models and the true cluster mass. The color code of different models in 
described in Sect. [4] Left panels refer to mass estimates performed using the entire sample of giants, right panels refer to estimates 
performed applying a radial selection to the sample. The epoch of core-collapse in simulation W5rhlR8.5 is marked by an arrow in 
bottom panels. 


neutron stars, binaries, a sub-system of stellar mass black 
holes or an intermediate-mass black hole, etc.). The results 
of the above test are described in the next subsections. 


4.1 Mass estimate 

In Fig. [5] the ratio between the mass estimated by the best- 
fit analytic models and the true mass of the simulation is 
plotted as a function of the cluster age. By analysing sim¬ 
ulation W5rhlR8.5 it appears that, after the first snapshot 
of the simulation, single-mass models systematically under¬ 
estimate the cluster mass up to 50% and improve their esti¬ 
mates at the end of the simulation. This is due to the quick 
relaxation reached by this simulation after few Myr: as al¬ 
ready explained in Sect, [l] in single-mass models all stars 
have a single kinetic temperature. When a significant level 
of relaxation is established massive stars (used to sample 
both the number density and the velocity dispersion of the 
cluster) have a smaller velocity dispersion with respect to 
other cluster stars, mimicking the behaviour of stars in the 
potential of a less massive cluster. Such an effect is less ev¬ 
ident when many mass bins are used to derive the number 
density profile, since in this last case the derived mass dis¬ 
tribution (and consequently the gravitational potential) is 


more similar to the true one. In the advanced stages of evo¬ 
lution the masses estimated by single-mass models increase 
and exceed the real value at t >10 Gyr. This trend occurs 
because of the flattening of the mass function due to the 
preferential escape of low-mass stars: in this stage the mean 
mass of stars increases and becoming more similar to that 
of the stars used as tracers of the cluster potential (see also 
SGlfff. Regarding multi-mass models with a = 1 they pro¬ 
vide a better estimate of the mass after ~2 Gyr, although 
they generally overestimate the cluster mass by 10-20%. A 
counterintuitive evidence is that models whose density pro¬ 
file is calculated only from the most massive bin perform 
slightly better than those where more bins have been con- 

5 A meaningful comparison between these two works can be made 
by comparing the M es t/Mtrue values of the single-mass King 
(1966) model fit made using only the most massive bin (magenta 
lines in Fig. 0 at ages >10 Gyr with the M 0 5 s /M plotted in the 
central panel of Fig. 3 in SG15 at a metallicity of [Fe/H]=-1 (note 
that their ’’Flat IMF” values must be divided by 0.41 to convert 
Mssp i n t° M). In both works masses appears to be underesti¬ 
mated by 10-25% (this work) and by 0-40% (SG15) with larger 
differences occurring when steeper MF are considered. Small dif¬ 
ferences are likely due to the different stages of dynamical evolu¬ 
tion of the clusters considered in the two works. 
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Figure 6. Ratio between the potential estimated by different 
models and the true potential of the snapshot at 11 Gyr of simula¬ 
tion W5rhlR8.5. The color codes of different models are described 
in Sect. [I] 

sidered. To understand the reason of such an occurrence we 
plot in Fig. [6] the ratio between the potential estimated by 
the various models and the true one in the snapshot at 11 
Gyr. It is apparent that while all models generally reproduce 
the potential at large distances from the cluster center, the 
central potential is always missed by models. This region is 
indeed particularly difficult to be reproduced because mas¬ 
sive remnants or binaries can sink into the very central part 
of the cluster producing a cusp-like shape not reproduced 
by modelfl We repeated the mass estimate by excluding 
the innermost region at r < rn/ 2. In this last case multi¬ 
mass models adequately reproduce the cluster mass during 
the entire evolution (within 10% except in the initial stages) 
regardless of the adopted range of masses used to compute 
the density profile. 

A different situation can be noted by analysing simu¬ 
lation W5rhll.5R8.5: here single-mass models always well 
reproduce the cluster mass, while multi-mass ones tend to 
overpredict the cluster mass by ~ 10%. Before 2 Gyr in¬ 
stead, multi-mass models assuming a = 1 significantly over¬ 
estimate (by more than 50%) the mass. This is because in 
these early stages of evolution of this simulation the cluster 
is still dynamically young and the low velocity dispersion 
of massive stars predicted by models is compensated by an 
increased estimated mass. The behaviour at t >5 Gyr is 
instead linked to the heating produced by the tidal field. 
Indeed, because of the strong tidal field experienced by the 
cluster in this simulation, many stars reach an energy level 
comparable to the potential at the tidal radius, but they take 

6 This effect is not removed by the adoption of a bin for massive 
remnants since it is produced by very few objects with mass larger 
than the adopted bin mean mass. 


some dynamical time to escape (Fukushige & Heggie 2000; 
Baumgardt 2001). During this phase they remain trapped 
within the cluster on energetic orbits at large radii with an 
increased velocity dispersion. In fact, by excluding the out¬ 
ermost region of the cluster (at r > rh where the orbits of 
these stars are preferentially confined) the estimated masses 
decrease and bring back the estimates of multi-mass models 
on the right value. It is worth noting that the best estimate 
is obtained constraining the value of a with two measures 
of the MF at different distances from the cluster center. 

4.2 Mass Function estimate 

In Fig.[7]the differences between the true and the estimated 
global MFs are shown for various assumptions and snapshots 
of the two considered simulations. Of course, only multi¬ 
mass models have been considered in this case. It is appar¬ 
ent that the various assumptions on the number of bins used 
to construct the density profile and the choice of the correct 
value of a do not produce significant differences in the de¬ 
termination of the MF, providing during all the evolution a 
good estimate of the MF slope (within A omf < 0.05). A 
systematic trend is instead visible in simulation W5rhlR8.5 
when the global MF is measured by correcting that esti¬ 
mated in the cluster center: in this case a slope shallower 
by A omf ~ 0.4 is estimated. This is a consequence of the 
insufficient level of relaxation of multi-mass models which 
is particularly evident in the cluster center (see Sect. IXT1) . 
Consequently, the lack of low-mass stars caused by mass 
segregation is erroneously interpreted as a real deficiency of 
stars and therefore lead to an underestimate of the MF slope. 
This bias is instead not visible in simulation W5rhll.5R8.5 
(where this last approach is anyway the most uncertain) 
because in this simulation multi-mass models better repro¬ 
duce the relaxation of the system. The difference between 
the actual global MF and that measured at the half-mass 
radius (without applying any correction) is also plotted in 
Fig. El for comparison. It is apparent that also in this case, 
no significant differences are noticeable. 


5 SIMULATION OF GAIA OBSERVATIONS 

In the previous sections we estimated the masses and the 
MFs of simulated clusters through the comparison with 
steady-state analytic models using as mass tracer the LOS 
velocity dispersion of the most massive cluster stars. In gen¬ 
eral, good estimates can be obtained with the desired level of 
accuracy provided that radial velocities for a large enough 
number of stars are available. The use of proper motions 
can improve the accuracy of the above estimates increas¬ 
ing the number of velocities used to constrain the cluster 
mass. However, the invaluable improvement provided by 
proper motions resides in the possibility of i) testing the level 
of relaxation by measuring the velocity dispersion of stars 
with different masses, and ii) measuring the characteristics 
of anisotropy. For the first task it is necessary to measure 
proper motions for MS stars down to a few magnitudes be¬ 
low the turn-off point. The faint limiting magnitude and the 
high angular resolution needed for this task can be achieved 
only with HST and will be soon available (see Bellini et al. 
2014). The second task can be instead obtained even using 
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Figure 7. Difference between the actual global MFs and those estimated using various assumptions for the location of the deep field 
range and the mass range of stars used to compute the density profile for the two considered simulations. The color code is described in 
Sect. [4] Red lines in both panels indicate the differences with respect to the MF observed at the half-mass radius without applying any 
correction. 


proper motions for only the brightest massive stars. Gaia 
will provide proper motions for stars brighter than g < 20 
over the full sky, sampling also the RGB stars of a number 
of nearby GCs. An, Evans & Deason (2012) showed that the 
proper motions provided by Gaia will allow an accurate es¬ 
timate of masses in GCs up to distances of ~20 kpc. In this 
section we want to test if a level of anisotropy like the one 
developed by the simulations presented here can be detected 
with the accuracy of Gaia data. 

In principle, we could use for this test the snapshots 
of the performed N-body simulations. However, the simula¬ 
tions considered here have masses ~10 times smaller than 
a typical GCs, so the velocity dispersions would be sig¬ 
nificantly smaller and the crowding conditions would be 
less severe than in a real GC. For this reason we simu¬ 
lated a set mock clusters from the potential and the phase- 
space distribution of the best-fit anisotropic multi-mass 
King-Michie model of the snapshot at 11 Gyr of the sim¬ 
ulation W5rhlR8.5, but assuming cluster masses between 
4.5 < log(M / Mq) < 6. Positions and velocities have been 
projected and transformed in angular distances and proper 
motions assuming various heliocentric distances. Gaia mag¬ 
nitudes have been calculated by interpolating the particle 
masses through an isochrone by Marigo et al. (2008) with 
suitable age and metallicity. A reddening of E(B — V) = 0.1 
and the extinction coefficients by Chen et al. (2014) have 
been assumed. Errors on proper motions have been added 
using the software PyGaia provided by the Gaia Project Sci¬ 
entist Support teanqj as a function of magnitude, color and 
assuming a K spectral type (suitable for RGB stars). These 
information are computed as all sky average, neglecting effi¬ 
ciency variations on the sky, and assuming the after-launch 
throughputs and error curves. Because of the angular size 
of the readout window, proper motions can be computed 
only for relatively isolated stars. Here we consider only stars 

7 http://github.com/agabrown/PyGaia 


with no neighbors brighter than 2.5 magnitudes with re¬ 
spect to the target g magnitude within 3.54" (see Pancino 
et al. 2013). Gaia will also provide radial velocities from 
low-resolution spectra. However, the estimated velocity er¬ 
rors and the limiting magnitudes are of lower quality with 
respect to those of available spectrographs mounted on 8m- 
class telescopes, so we considered errors of 0.5 km/s on radial 
velocities only for stars brighter than g < 17 assuming they 
come from a different source. 

The anisotropy parameter (3 is a function of the ratio 
between the tangential and the radial components of the 3D 
velocity dispersion (a 2 and a 2 , respectively). These quanti¬ 
ties are related to the velocities along the LOS ( vlos ) and 
projected on the plane of the sky ( v acos s and vs ) by Eq. [5] 
where 


x = (a — ao) cosS 
y = 5 - 5 0 

R — \J x 2 - + y 2 

r = yjx 2 T y 2 T z 2 
xVacosS T yrs 



Qo and <5o are the right ascension and declination of the 
cluster center and 2 is the distance along the LOS from the 
plane passing through the cluster center. As the distance 
along the LOS is unknown it is not possible to determine /3 
directly from observations. On the other hand, it is possible 
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Figure 8. Significance ((/3)/cr^) of the anisotropy detection as a 
function of the cluster distance. Black, red, green and blue lines 
indicate simulations of GCs with masses log(M/M q)= 4.5, 5, 5.5 
and 6, respectively. Solid lines indicate estimates done using only 
stars with all the three velocity components, dashed lines indicate 
estimates done using all stars with available proper motions. 


to combine Eq.s[4]to obtain 


P{R) = 


L 


rt n(r)<j^./3r 
R y/r 2 -R? 


dr 


L 


rt ra(r)crpr 
R 


dr 


= 1 - 


a r' + a LOS 


Since the maximum density and radial velocity dispersion 
in the above equation occurs at (r = R, z = 0), the above 
quantity is close to ft(R). 

One thousand synthetic observations have been simu¬ 
lated for each choice of the cluster mass and distance and the 
mean value of {/?') and its standard deviation api have been 
calculated for those stars with both radial velocities and 
proper motions estimates. The /3 1 value of each synthetic ob¬ 
servation has been calculated using the maximum-likelihood 
algorithm described by Eq. 0 to calculate ctlos , o>/ and 
ay where individual errors on velocities (ei) have been 
added to the intrinsic velocity dispersion ( (Jintr ) assuming 
of = crjntr + e? ■ The significance level (defined as (/ 3')/api ) 
is plotted in Fig. [8] as a function of the adopted cluster dis¬ 
tance for different simulated masses. It can be seen that for 
low/intermediate-mass clusters (log(M / Mq) < 5.5) it is not 
possible to detect with a good level of significance (> 3) the 
anisotropy, while a better signal is detectable for massive 
(log(M /Mg) ~ 6) clusters up to distances of ~15 kpc. 

In the above test we used only stars brighter than g < 17 
with both proper motions and LOS velocities. These two 
observables are however derived from completely different 
techniques and are subject to different sources of uncertain¬ 


ties. For instance, proper motions are measured in terms 
of angular motions and require an estimate of the cluster 
distance to be converted into physical units, at odds with 
radial velocities coming from the observed Doppler-shift of 
spectroscopic lines. Moreover, because of the different tech¬ 
niques, the estimated velocity errors (essential to properly 
calculate intrinsic velocity dispersions) can present different 
systematic biases. So, we considered also the case where only 
proper motions are available. In this case we considered the 
quantity 



Although the part of the information coming from radial 
velocities is lost in the above parameter (/?" is a worse esti¬ 
mator of /3 than /3 1 ), its corresponding uncertainty is smaller 
since i) it is given by the propagation of a smaller number 
of uncertainties, and ii) many more stars can be used to 
compute P" (down to g < 20). The significance level of this 
approach for the different adopted masses and distances is 
overplotted in Fig. [8] It is apparent that this last approach 
provides a larger significance level: although at large helio¬ 
centric distances significant detections remain restricted to 
the most massive GCs, at distances <6 kpc it seems pos¬ 
sible to detect signatures of anisotropy also in less massive 
clusters (log(M /Mq) > 5). 

It has been shown (see Sect. 0 that the level of 
anisotropy changes with radius both in models and in N- 
body simulations. So, the values of ft' and /3”, which are av¬ 
erages over stars sampled at different distance from the clus¬ 
ter center, cannot provide a fair indication of the character¬ 
istics (strength, radial profile, etc.) of the cluster anisotropy. 
A lager number of proper motions covering the entire cluster 
extent are needed for this purpose. 


6 SUMMARY 

The comparisons presented in Sect. 0 show that multi-mass 
models can qualitatively reproduce the radial distribution 
and velocity dispersion of stars with different masses in GCs 
subject to frequent collisions during a wide range of their dy¬ 
namical evolution. In particular, the formulation by Gunn & 
Griffin (1979; a=l) is able to provide a fair representation 
of N-body simulations after a timescale comparable to the 
half-mass relaxation time. By looking at Milky Way GCs 
this condition is satisfied by all but three GCs (namely uj 
Gen, Pal 14 and NGC2419; see McLaughlin & van der Marel 
2005; Marin-Franch et al. 2009). During the initial half-mass 
relaxation time multi-mass models overestimate the status 
of relaxation predicting an unrealistic level of mass segrega¬ 
tion. In this case, a simple generalization of the Gunn & Grif¬ 
fin (1979) formalism (adopting At oc m“) provides models 
accounting for intermediate levels of relaxation which better 
reproduce the simulation properties. As a note of caution, 
consider that the N-body simulations used here start with¬ 
out any degree of primordial mass segregation which appears 
to be present in many young massive clusters (de Grijs et al. 
2002; Frank, Grebel & Kiipper 2014). After many relaxation 
times the analysed simulations undergo strong mass segre¬ 
gation which is not reached by multi-mass models which 
slightly underpredict the mass segregation present in the 
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cluster. Both the considered simulations quickly develop two 
opposite velocity anisotropy trends which cannot be repro¬ 
duced by analytic models. In the Gunn & Griffin (1979) for¬ 
malism only radially anisotropic velocity distributions can 
be simulated in King-Michie models, but alternative formu¬ 
lations exist to account for tangentially anisotropic ones (see 
e.g. Weinberg 1991; An & Evans 2006). However, the simple 
dependence of the King-Michie distribution function from 
angular momentum per unit mass cannot reproduce the wide 
range of trends of the anisotropy parameter /3 with stellar 
masses observed in N-body simulations. New models need 
to be developed to model this quantity in real GCs. 

Besides the above considerations, we tested the accu¬ 
racy of simple fits of the generally available observables (pro¬ 
jected number density profiles from RGB+bright MS stars, 
LOS velocities of RGB stars, MF estimated in a spatially re- 
striced region) with analytic models in the estimate of mass 
and MF. It appears that the adoption of single-mass mod¬ 
els to model GCs after many half-mass relaxation times (as 
widely done in many past studies because of the lack of in¬ 
formation on the cluster MF) can lead to underestimates of 
their masses up to ~50%. This is likely the reason why dy¬ 
namical masses estimated in early studies resulted smaller 
than luminous ones (Mieske et al. 2008). The same con¬ 
clusion has been reached by SG15 who quantified the bias 
in the mass estimated from integrated properties of a set 
of multi-mass GCs models if mass segregation is neglected 
and found that this effect can lead to a mass underesti¬ 
mate ranging from 0.25 < M 0 bs/M tr ue < 1 depending on 
the metallicity, the MF and the retention fraction of dark 
remnants (see their Fig. 3). In the same way, the use of 
multi-mass King-Michie models (with a = 1) in dynami¬ 
cally young GCs produces the opposite bias leading to a 
significant overestimate of the cluster mass. In this case, the 
adoption of models accounting for an intermediate level of 
relaxation (0 < a < 1) provides better results, although it 
requires the measure of the MF in two different radial ranges. 
Both the central cluster region and the outskirts are regions 
sensitive to the uncertain distribution of heavy objects (like 
neutron stars, black holes and binaries) and tidal heating. In 
this respect, tidal heating can be effective in increasing the 
velocity dispersion even close to the half-mass radius and 
therefore the estimated mass (as already found by Kiipper 
et al. 2010). An analysis devoted to the derivation of GCs 
masses should therefore be restricted to a relatively narrow 
portion of the cluster (rn/2 < r < rh) to minimize the 
impact of the above mentioned effects. In the simulations 
considered here, the adoption of different samples of stars 
and the location of the deep field used to measure the MF 
do not significantly affect the mass estimation if multi-mass 
models are considered. This is due to the relatively good 
representation of the mass segregation process provided by 
multi-mass models. 

The estimate of the MF is instead less affected by the 
uncertainties in the density profile and the level of relax¬ 
ation. However, in strongly mass segregated clusters the MF 
slope can be underestimated by Aojuf ~ 0.4 if the local MF 
is measured at the cluster center. This is due to the inade¬ 
quacy of multi-mass King-Michie models in reproducing the 
relative distribution of low/high mass stars in the central 
region of such relaxed GCs. Note that both the most ex¬ 
tensive works on the estimation of the MF slopes (Paust et 


al. 2010; Sollima et al. 2012) are based on MF estimated 
from the same dataset of HST observations in the central 
region of GCs being potentially affected by this bias. This 
effect is however negligible in those dynamically young clus¬ 
ters where multi-mass King-Michie models do a good job in 
reproducing the radial distribution of different masses. It is 
not clear how many GCs of the above mentioned studies can 
suffer from this bias. On the other hand, the MF measured 
at the projected half-mass radius (without the application of 
any correction) is already a fairly good representation of the 
cluster global MF (see also Pryor et al. 1986, Vesperini & 
Heggie 1997; Baumgardt & Makino 2003). So, studies based 
on MF estimated at the half-mass radius (like those by Pi- 
otto et al. 1997 and Piotto & Zoccali 1999) can provide good 
estimates of the global MF, provided that a good estimate 
of the half-mass radius is available. 

A limitation of the approach adopted here is that the 
clusters simulated in our N-body simulations are a factor of 
5-10 smaller than real GCs. Since the effects of stellar evolu¬ 
tion, relaxation and tidal effects on the dynamical evolution 
of the cluster act on different timescales and their efficien¬ 
cies depend on the cluster mass in different ways, it is not 
possible to adopt some sort of scaling to GC-size objects 
(see Baumgardt 2001). We recall that the simulations anal¬ 
ysed here consider two extreme situations: while simulation 
W5rhlR8.5 is a system undergoing quick core-collapse, sim¬ 
ulation W5rhll.5R8.5 is dynamically young and subject to 
a relatively strong tidal field during all its evolution. Thus, 
we expect that these simulations bracket the whole sample 
of possible evolutions occurring in GCs. Anyway, any com¬ 
parison with N-body simulations cannot be used to draw 
general conclusions on the dynamical properties of real GCs. 
Further studies on real GCs stars are needed to assess the 
validity of the results presented here. 

We explored the possibility to detect the presence and 
the characteristics of anisotropy in GCs using the proper 
motions provided by Gaia. We find that moderate levels of 
anisotropy (like those naturally developed by N-body sim¬ 
ulations) can be detected in massive ( log(M/MQ ) ~ 6) 
GCs at distances <20 kpc and in intermediate-mass GCs 
{log(M/M q) > 5) GCs at d < 6 kpc. By looking at the Har¬ 
ris catalog (Harris 1996; 2010 edition) there are 25 GCs in 
this radial/mass range. Unfortunately, with the information 
provided by Gaia only, it will not be possible to sample the 
number of objects needed to construct anisotropy profiles 
even for the nearest GCs. Indeed, the limiting magnitude 
and the spatial resolution of Gaia will allow to obtain proper 
motions with errors better than <j v < 2 km/s only for RGB 
stars in the outskirts of GCs. Nevertheless, Gaia data will 
be complemented by HST proper motions of stars in the 
central regions of GCs thus providing the needed coverage 
for these kind of studies. 
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